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Abstract The Moran model with recombination is considered, which describes the evolution of the genetic 
composition of a population under recombination and resampling. There are n sites (or loci), a finite 
number of letters (or alleles) at every site, and we do not make any scaling assumptions. In particular, 
we do not assume a diffusion limit. We consider the following marginal ancestral recombination process. 
Let S = {1,..., n} and A = {Ai ,..., Am} be a partition of S. We concentrate on the joint probability 
of the letters at the sites in Ai in individual 1, ..., at the sites in Am in individual m, where the 
individuals are sampled from the current population without replacement. Following the ancestry of 
these sites backwards in time yields a process on the set of partitions of S, which, in the diffusion limit, 
turns into a marginalised version of the n-locus ancestral recombination graph. With the help of an 
inclusion-exclusion principle, we show that the type distribution corresponding to a given partition may 
be represented in a systematic way, with the help of so-called recombinators and sampling functions. The 
same is true of correlation functions (known as linkage disequiiibria in genetics) of all orders. 

We prove that the partitioning process (backward in time) is dual to the Moran population process 
(forward in tim e), where the sampling f unction plays the role of the duality function. This sheds new light 


on the work of Bobrowski et al. ( 201011 . The result also leads to a closed system of ordinary differential 


equations for the expectations of the sampling functions, which can be translated into expected type 
distributions and expected linkage disequiiibria. 

K©ywords: Moran model with recombination, ancestral recombination process, linkage disequiiibria, 
Mobius inversion, duality 

MSC 2010 Subject classification: 92D10, 60J28. 


1 Introduction 


Models that describe the evolution of hnite populations under recombination are among the 
major challenges in population genetics. This article is devoted to the Mor an model with 
recombination (in con tinuous time), which is briefly described as follows (see Durrett 20081 : 


Bobrowski et al. 2O10l ). A chromosome is identified with a linear arrangement (or sequence) 


of n discrete positions called sites, which are collected in the set 5' = {l,...,n}. A site may 
be understood as a nucleotide site or a gene locus. We will throughout consider chromosomes 
as (haploid) individuals, that is, we think at the level of gametes (rather than that of diploid 
individuals that carry two copies of the genetic information). Site i is occupied by letter Xi G 
where Xj is a hnite set, 1 ^ i ^ n. If sites are nucleotide sites, a natural choice for each 
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Xj is the nucleotide alphabet {A, G,C,T}; if sites are gene loci, Xj is the set of alleles that 
can occur at locus i. The genetic type of each individual is thus described by the sequence 
X = {xi,X 2 , ■ ■ ■, Xn) E X;^ X • • ■ X X^ =: X, where X is the type space. Recombination means that 
a new individual is formed as a ‘mixture’ of an (ordered) pair of parents, say x and y. We will 
restrict ourselves to single-crossover recombination, that is, the offspring inherits the leading 
segment (up to site i, for some 1 ^ i < n) from the hrst and the trailing segment (after site 
i) from the second parent. The recombined type thus is (x<jj, y>i) ;= (xi,..., x*, y^+i,..., Vn)', 
we say that a crossover has happened between sites i and i + 1. The sites that come from the 
paternal and the maternal sequence, respectively, dehne a partition ^ of 5 into two parts (we 
need not keep track of which part was ‘maternal’ and which was ‘paternal’). All partitions of S 
into two ordered (or contiguous) parts {A = {{1, 2 ,... ,i}, {i + 1 ,..., re}}, i € S \ {re}) can be 
realised, via a single crossover event. Altogether, whenever an offspring is created, its sites are 
partitioned between parents according to A with probability r^, where ^ 0, X]^e 02 (S) Gt ^ 
and 02 ( 8 ) is the set of all ordered partitions of S into two parts. Let us note that, due to the 
one-to-one correspondence between elements of S' \ {re} and those of ©2(5), the specification of 
the simply means that a crossover probability is associated with each site in S\{re}. The sum 
S, 4 g 02 (S) La probability that some recombination event takes place during reproduction. 

With probability = 1 — X^yig 02 (S) LA’ there is no recombination, in which case the offspring 
is the full copy of a single parent. We write 0^2(S) := ©2(5) U {S} for the set of ordered 
partitions in to at most tw o parts. The collection {r 4 }_ 4 go< 2 (S) known as the recombination 
distribution (Burger 20001 . p. 55). 


Consider now a population of a constant number of N haploid individuals (that is, gametes), 
which evolves as follows (see Figure [T]). Each individual has an exponential lifespan with 
parameter 1 (this choice of the parameter is without loss of generality; it simply sets the time 
scale). When an individual dies, it is replaced by a new one as follows. First draw a partition A 
according to the recombination distribution. Then draw |.4.| parents from the population (the 
parents may include the individual that is about to die), uniformly and with replacement, where 
|.4.| is the number of parts in A. If \A\ = 2, the offspring inherits the leading segment of A 
from the hrst and the trailing segment from the second parent, as described above. If |^| = 1 
(and thus A = {S}), the offspring is a full copy of a single parent (again chosen uniformly from 
among all individuals); this is called a (pure) resampling event. All events are independent of 
each other. 


Note that it may seem biologically more realistic to draw two parents without replacement. 
However, assuming sampling with replacement entails signihcant simplihcations, and yields 
the same process as sampling without replacement with a slight change in the recombination 
distribution. More precisely, since drawing the same individual twice means that the offspring 
is a full copy of this single parent, our process agrees (in distribution) with the analogous 
process without replacement if is replaced by r_^{N — 1)/N for all A E 02(*S') (and r^^^ is set 
accordingly). 


The model will be described more formally later. For now, let us summarise the two main 
lines of research in this context. On the one hand, there has been considerable interest in how 
the composition of the population evolves over time, and, in particular, how the correlations 
betwe en sites (kno wn as linkage disequilibria) will develop; see the overviews in i Hein et al 


( 2005 . Chap. 5.4), Durrett ( 20081 . Chap. 3.3 and 8.2), or Wakelev ( 20091 . Chap. 7.2.4). Since 
there is no mutation, a single type will go to fixation in the long run, that is, the entire 
population will ultimately consist of this single type. In the absence of recombination, this 
will be one of the types initially present, and it is well known that the fixation probability for a 
given type equals its initial frequency. If there is recombination, the type that ultimately wins 
can also be a newly-composed type, but little is known about the fixation probabilities of the 
many possible types. The explicit development over time is even more challenging, due to an 


- 2 - 
























T I 


Figure 1: Snapshot of a Moran model realisation with N = 5 individuals. For example, in the first event, 
individual 3 dies and is replaced by a recombined copy of individuals 2 and 3. The last line shows the composition 
of the population at the final timepoint, T. 


intri cate interplay of resampUng and recombination. It is usually approached fo rwards in time 


e.g., lOhta and Kimural (Il969li. 
Chap. 8.2), or lBobrowski et al 


Song and Song (2007), Baake and Hermi ( 2008l h Durrett ( 20081 . 
( 2O10l ). In the deterministic limit, which emerges when N 


without rescaling of the or of time, the population is described by a system of ordinary 
differential equations, again forwards in time. This system has an explicit solution, both for 
the type distribution and for correlat ion functions of all orders, for an arbitrary number of sites 
( Baake and Baake 20031: Baake 2005 ) . Thi s also provides a decent approximation for large but 
finite populations ( Baake and Her ms 200a), but dealing appropriately with the stochasticity of 
finite populations remains a major challenge. 

The second line of research is concerned with genealogical aspects and sampling formulae. Here, 
one starts from a sample taken from the present population and traces back the ancestry of the 
various segments the individuals are composed of. A major challenge lies in the calculation of the 
probabilities for the type dist ribution of a random sample, that is, one aims at so-called sampling 
formulae, see [Durrett] (j2008l . Chap. 3.6). These questions are naturally approached backwards 
in time. Usually, one employs the diffusion (or weak recombination) limit, that is, time is sped 
up by a factor of N, followed by ^ oo such that Nrj^^ Qj^, A ^ 02(5'). Obtaining sampling 
formulae is tied to the situation in which the population has reached a stationary state; even 
this case is very hard to treat, and coping with time dependence seems to be hopeless. 


The aim of this article is to build a bridge between these two lines of research. We will explore 
the type distribution an d the correlations over time, in the stocha stic setting. A starting point 
will be a recent paper by Bobrowski. Woidvia. and Kimmel ( 2O10l ). who approach this question. 
Their setting is entirely forwards in time, which effectively hides some of the underlying structure. 
In contrast, we will proceed backwards in time and provide a genealogocial approach for the 
analysis of correlations. The crucial notion in this context will be that of duality between the 
original Moran model forward in time and a suitable ancestral process that follows back the 
ancestry of selected segments fr om today’s population. This will also shed new light on the 


results of Bobrowski et al. ( 2nifll h In order to keep the approach as general as possible, we will 


throughout adhere to the original (finite N) model, without taking any limit, but will discuss 
the various scalings and limits where appropriate. 


The paper is organised as follows. In Section [2j we start by collecting some important facts about 
partitions and Mobius functions. We then introduce the model more formally and motivate 
our genealogical approach in Section [3l which may be considered a marginal version of the 
usual ancestral process with recombination. In Section HI we describe our ancestral process, 
which is a partitioning process that keeps track of how the ancestral material is partitioned 
between individuals. In Section HI we introduce a systematic description via recombinators, 
which describe the action of recombination on a population and have proved very useful in the 
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deterministic setting. We complement them here by sampling functions, which are additionally 
required for hnite populations. The collection of sampling functions will be crucial since it 
will also serve as duality function in Section [6l where the duality between the Moran model 
forward in time and the partitioning process backward in time is proved. This proof, at the 
same time, yields a differential equation system for the expectations of the sampling functions, 
which are the building blocks for the linkage disequilibria. In Section [71 we apply our results to 
the cases of two and three sites. We will see that the expected linkage disequilibria (of second 
and third order) decay exponentially even in the presence of resampling, and identify further 
linear combinations of expected sampling functions that decay exponentially. For two sites, we 
also obtain the explicit time course for the expected composition of the population, and, at the 
same time, the hxation probabilities of the various types. 


2 Preliminaries: Partitions and Mobius functions 


Working with partitions will be essential to our approach, and we will rely throughout on the 
powerful concept of Mobius functions and Mobius inversion. Let us briefly collect the basic 
notions and standa r d resu lts; more b ackground ma terial as well as th e proofs may b e found in 
Rota (1964), Berge ( 1971 . Chap. 3.2), Aigner ( 1979 . Chap. I,II,IV) and Stanley ( 19861 . Chap. 3). 


Partitions. Let IT be a finite, nonempty, totally ordered set, such as a finite subset of N; later, 
W will be 5 or a subset thereof. Let P = P(IT) be the set of partitions of W. We write such a 
partition as A = {Ai,, A^}, where Aj A ^ for fol f and Aj n = 0 for all j ^ k together 
with Ai U • • • U Am = W. We call Aj a block (or part) of A and m = |A.| is the number of blocks 
in A. 

We say that a partition A = {Ai,..., Am} of W is ordered (or contiguous, or an interval 
partition) if every Aj is ordered in W, that is, Aj = {x ^ W \ minAj ^ x ^ maxAj}. For 
example, if VF = {1,2,5, 7,9}, then {{1, 2, 5}, {7, 9}} is ordered, but {{1, 2, 7}, (5,9}} is not. 
The set of all ordered partitions of W is denoted by 0(IT), the set of all ordered partitions of 
W into (exactly) two parts is 02(IT), and the set of all ordered partitions of W into at most 
two parts is 0^2(IF)- 

For a given partition A = {Ai,..., Am} of W, let M := (1,2,..., m} = M{A) and, for J C M, 
we define Aj := {Aj}j^j and Aj := Uj^jAj. Aj is a partition of Aj. In particular. Am = A, 
Am = IF) ~^{j} = AM\{j} = -4. \ {Aj}, for any j G M. Note that M depends on 

A, but we suppress this dependence when there is no risk of confusion. We will throughout 
abbreviate J \j := J \ {j} and J U A: := J U {k}. 

The natural ordering relation on P(IT) is denoted by where A ^ B means that A is a 
refinement of B, that is, every block of A is a subset of a block of B; equivalently, B is a 
coarsening of A. A B means that A ^ B and A A Together with the resulting partial 
order, P(IT) is a poset and, in particular, a finite lattice. P(IF) has a unique minimal or finest 
partition, which is denoted as 0 = {{x} | x € IF}; likewise, there is a unique maximal or coarsest 
one, namely 1 = {IF}. 

When U and V are disjoint (finite) sets, two partitions A G P(t/) and B G P(F) can be joined 
into A U B to form an element of P(t/UF). Furthermore, if 17 C IF, a partition A G P(IF), 
with A = {Ai,..., Am} say, defines a unique partition of U by restriction. The latter is denoted 
by A\u, and its parts are precisely all non-empty sets of the form Aj n 17 with 1 A i A mn. In 
particular, 1|^ is the coarsest element in P(17). For two partitions A and B, the least upper 
bound will be denoted by A V B, namely the hnest partition C for which A ^ C and B ^ C. 
Analogously define the greatest lower bound of A and B hy A A B. 
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Example 1 . Consider the two partitions A = {{1,3,4}, {2, 5}} and B = {{1,4}, {2,3}, {5}} of 
W = {1,..., 5} together with a subset U = {1, 2,4} of W. Then AaB = {{1,4}, {2}, {3}, {5}}, 
AyB = {{1,... , 5}}, and A\u = {{1,4}, {2}}. 0 


Mobius functions on the poset of partitions and Mobius inversion. The Mobius function 
of a poset is a general and powerful tool in discrete mathematics. It may be considered as 
a systematic way of implementing the inclusion-exclusion principle. We rely on it in two 
contexts in this article: First, we use it to turn sampling without replacement into sampling 
with replacement, and vice versa. Second, we need it to turn type frequencies into linkage 
disequilibra. 


For b ackground mater ial, we r e fer th e reader to iRotal (|l964l l. iBergd (jl97ll . Chap. 3.21. lAigner 
( 1979 . Chap. IV) and Stanley ( 198fil . Chap. 3). Let us only summarise here that the Mobius 
function p, is defined inductively for all A ^ B € P(VF) via 


fi{A,A) = l and fj.{A,B) = — ^ fi{A,C), (1) 

A^C<B 


where the underdot indicates the summation variable. Let A ^ B a P(1T), with m = \B\ the 
number of blocks in B, and Uj the number of blocks of A within block Bj of B, that is, nj is the 
number of blocks in A\^., 1 ^ j ^ m. The Mobius function of the pair {A,B) is then given by 

m m 

pi{AB) = n KAb,A\b,) = - 1)!. 

1=1 1=1 

The following identity for sums of Mobius functions 

^ = ( 2 ) 

1°’ Otherwise, 

follows immediately from the definition (Jll and will be used frequently in what follows. 

We can now state the fundamental Mobius inversion principle. Let / and g be mappings from 
P(1T) to C which are, for all A G P(1T), related via 

g{A) = ^ fiB). (3) 

B^A 

Then, this can be solved for / via the inversion formula 

f{A)=^f,{A,B)giB). (4) 

B)^A 

More precisely, this is inversion from above. An analogous formula applies for inversion from 
below; this relies on refinements rather than coarsenings, with replaced by in ([2D and (HD. 
It is important to note that Mdbius inversion is not restricted to functions; it also applies to 
bounded operators. 


3 The model and the genealogical approach 

In this section, we define the model formally and motivate our genealogical approach. 
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3.1 The Moran model with single-crossover recombination 


We identify the population at time t by a (random) counting measure on X. Namely, Zt{{x}) 
denotes the number of individuals of type x € X at time t, and Zt{A) := Xt({x}) for 

^ C X; we abbreviate Zt{{x}) as Zt{x). If we define 6 x as the point measure on x (i.e., 
^x{y) = Sx,y for x,y G X), we can also write Zt = J2xex ^t{x)6x- Since our Moran population 
has constant size N, we have ||Xt|| = N for all times, where \\Zt\\ := ~ is the 

norm (or total variation) of Zt- 

So, {Zt}t^o is a Markov process in continuous time with values in 

^:= {ze {0,...,iV}W I ||z|| =iV}, (5) 


where |X| is the number of elements in X. We will describe the action of recombination on 
(posit ive) meas ures with the help o f so-c alled recombinators as introduced bv iBaake and Baake 
( 2003l i: see also Baake and Hernii (2008) for a pedestrian introduction. Let A1+(X) be the set 
of all positive, finite measures on X and we understand W(+(X) to include the zero measure. 
Define the canonical projection vr^-: X —)■ =• ^i-, for I ^ S,hy 'Xj{x) = {xi)i^i as usual. 

For uj G A1_|_(X), the shorthand vr/.cu := oj o indicates the marginal measure with respect to 


the sites in / C S', where is the preimage of tt/. The operation . (where the dot is on the 
line and should not be confused with a multiplication sign) is known as the pushforward of u 
w.r.t. TTj. In terms of coordinates, the definition may be spelled out as 


[tTj.Uj) ( xj ) = 00 O TTj^ ( xj ) = 0 o[{x G X I TTj{x) = Xj}) , Xj G X/. 


Note that tt^.lo = ||a;|| and vr^.o; = to. 

Consider now A = {{1, 2,... , i}, {i + 1,..., n}} G ©2(5) and oo G WI+(X) \ 0, and define the 
projective reeombinator as 

^ ® ( 6 ) 

Moreover, we set R^{oj) := a;/||w||. R^{oo) is a probability measure for all oo G WI+(X)\0, where 
the zero measure is excluded to make it well-defined. In words, R^ turns oo into the (normalised) 
tensor product (or product measure) of its marginals with respect to the blocks in A. Writing 
out ([U]) in terms of coordinates gives 

{R^{UJ)){X) = 

D,) K 


where means marginalisation. If w = 2 ; is the current population, then R^(z) is the type 
distribution that results when a new individual is created by drawing a leading and (possibly) a 
trailing segment (as encoded by X G 0<g2(<S')) from the current population, uniformly and with 
replacement. 


Remark 1. R^ is a projective version of the reeombinator defined by Baake and Baake ( 200,81 1 : 
it differs from the latter by a factor of I/||a;||. Clearly, both versions agree on the set of probability 
measures. As we shall see, the projective version is more suitable in the stochastic setting, while 
the original recombinators are better adapted to the deterministic situation. Since recombinators 
will only appear in the projective version in this article, we will drop the superscript and the 
speciheation ‘projective’ and call := R^ a reeombinator by slight abuse of language. <0 
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In Section [5l we will generalise the recombinators and learn more about their probabilistic 
meaning and mathematical properties. For the moment, let us use them to reformulate the 
Moran model with recombination in a compact way. Namely, since all individuals die at rate 1, 
the population loses type-7/ individuals at rate Zt{y). Each loss is replaced by a new individual, 
which is sampled uniformly from Rj^{Zt) with probability r^, A € 0 ^2{S). Therefore, when 
Zt = z, the transition to z + 5x — Sy occurs with rate 

X{z-y,x):= r_A{R^{z)){x)z{y). (7) 

.4e0^2(S') 


The summand for A = 1 corresponds to pure resampling, whereas all other summands include 
recombination. Note that A includes ‘silent transitions’ {x = y). 

Remark 2. We would like to mention that the model may alternatively be formulated in 
terms of reproducing individuals rather than dying individuals, as follows. Each individual 
reproduces at rate 1 and picks an A ^ 0^2(5') according to the recombination distribution. If 
A O 2 , the reproducing individual contributes the sites in one of the blocks in A and picks 
a random partner that contributes the sites in the other block to the offspring. If ^ = 1, 
the reproducing individuals contributes all sites. The offspring pieced together in this way 
replaces a uniformly chosen individual from the population. In this formulation, which is closer 
in spirit to the deterministic single-crossover model, offspring of type x are created at rate 
Nrj^[Rj^{Zt)){x) and replace an individual of type y with probability Zt{y)/N. This explains the 
different normalisation of the original recombinator, w hereas the additional factor of N = \\Zt\\ 
is absorbed in its definition in Baake and Baake ( 2003l l . The resulting transition rates, however, 
are again those in Eq. dZD- Here, we stay with the formulation that led to Eq. ([7|) in the first 
place, since it seems more natural for finite populations. <C> 


Let us summarise our model as follows: 


Definition 1 (Moran model with single crossovers). The Moran model with single crossovers 
is the Markov chain in continuous time {Zt}t^Q with state space E of Eq. Q and generator 
matrix A with nondiagonal elements 

A{z,z + w)= ^ X{z-,y,x), wAO, 

Sx —Sy =W 

for z € E, w ^ E — z (where E — z := {v \ z + v € E}) and A{z, z) = — ^ 71(2:, z + v). 

v^E—z: 

v^O 


Limits of the forward modei. Consider now the family of processes {zj:^'^}t^Q, N = 1,2,, 
where we add the upper index to indicate dependence on population size. Also consider the 
normalised version {Z^^^/N}t^Q-, zj:^'^/N is a random probability measure on X. Eor N ^ 00 
and without any rescaling of the or of time, the sequence {zj:^'^}t^Q converges to the solution 
of the deterministic single-crossover equation 

iot= ^ rj^{R_^{uJt) - ujt) ( 8 ) 

>le02(S) 


with initial value Wg, Wg a probability m easure, and we assume t hat limjv^cxD Xg'^^/E = Wg. 
This is a law of large numbers and due to lEthier and Kurta (jl985l . Thm. 11 .2.1). The precise 


statem ent as well as the proof are perfectly analogous to Proposition 1 in iBaake and Her ms 


(j2008l L which assumes a slightly different sampling scheme for recombination. We therefore 
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leave out the deta i ls her e. T he dete r minis tic single-crossover equation ([8|) was investigated by 
Baake and Baake (20^) and Baake ( 20051 1 . For comparison, note that, in view of Remark [21 
the probability in ([8]) is multiplied by the unit rate at which each individual reproduces, and 
this way turns into a recombination rate. 

The Moran model with recombination also has a well-known diffusion limit, which emerges when 
N ^ oo under Nr^ —>■ A € ©2(5), afte r a speedup of time by a factor of N. In the case 

of two loci and two alleles, this goes back to Ohta and Kimnra ( 1969l l: see also Dnrrett (20^, 
Chap. 8.2) for a modern e xposit ion. Two loci with an arbitrary (but finite) number of alleles are 
treated by Jenkins et al. ( 20151 1. This can be readily generalised to the case of a finite number 
of loci with a finite number of alleles, but we do not spell it out here, since we will not draw on 
the diffusion limit of the forward process later. 


3.2 The ancestral recombination process (ARP) and its marginal version 


In line with standard population-genetic thinking, we employ a genealogical approach by tracing 
back the ancestry of (parts of) the genetic material from a population at present that evolved 
according to the Moran model with single-crossover recombination. The standard genealogical 
approach for models with recombination is the ancestral recombination graph (ARC) first 
formulated bv lHudsonl (j 19831 ). Today, many different notions of ‘ARC’ are in use. We stick to the 
usual convention here that the ARG assumes the diffusion limit. Hudson’s o riginal version was 


for t wo loci, but multilocus generalisations ([Griffiths and Mar 


201 2h and continuous 
immediate. 


sequence versions (n —>■ oo, see, e.g., iDnrrettl I2nn8l . Ghao. 3.4) 


pram 


19961: Bhaskar and Song 


are 



Figure 2: A realisation of the full ancestral recombination 
process, starting from m = 3 individuals; ancestral material 
is shaded, non-ancestral material is indicated by thin 
horizontal lines. The mixed recombination-coalescence 
event indicated by dashed lines can only appear in the finite 
population recombination process (ARP). In the diffusion 
limit, and thus in the ARG, recombination and coalescence 
act in isolation. 


The ARG starts from a sample of individuals from the present population and follows their 
ancestry backwards. When a sequence (or a part of a sequence) experiences a recombination 
event, it branches into a leading and a trailing segment; when two (p arts of) sequences go back to 
a comm o n anc estor, there is a c oalescenc e event . For overviews see lHein et al.l (|2005l . Ghap. 5), 
Durrett (20^, Ghap. 3.4), or IWakelevI (|2nn9l . Ghap. 7.2). Mutation can be independently 


superimposed on the ARG, but will not be considered in this article. One is then interested in 
the full information on the sample, namely, the probabilities for all possible type distributions 
of the sample. The stationary state of the ARG may be characterised by a collection of 
so-called sampling recursions'., they may be solved analytically for tiny s ample s (le ading to 


explic it sampling formulae), or numerically for larger ones, see Golding ( 19841 ) . or Durrett 
([20081 . Ghap. 3.6). But feasibility is limited due to the enormous state space, even for small 
samples. Alternatively, one resorts to computationally intensive Monte-Garlo or importance 
sampling methods to simula te the ARG ( Griffiths and Marjoram 19961 : Wang and Rannala 20081 : 


Jenkins and Griffithsll201ll ). Recently, Song and coworkers discovered structural properties of 


the ARG that allow for an efficient combination of analytical and simulation techniques in the 


















































































































Figure 3: The marginalised version corresponding to the ARP in Figure [2] in which we only follow blocks of the 
partition (shaded), that is, block Ai is sampled in individual number i, 1 ^ i ^ m. Material that is ancestral to 
the sampled individuals, but not to the blocks considered, is shown as open rectangles (left). But since this is 
not traced back, it can be treated in the same way as material non-ancestral to the sampled individuals (right). 
Consequently, the sample will finally consist of the blocks of the partition only. 


regime of strong recombination ( Jenkins and Song 201(11 ): more precisely, they work in terms of 
expansions in 1 /g g ^ cx), where all gj^, with A € 02(5), are assumed to scale linearly with 
the common factor g. 

In contrast, we will work in the setting of both finite n and finite N. The corresponding ancestral 
recombination process (ARP), which is illustrated in Figure [3 is a finite-population version of 
the multilocus ARG. We then simplify matters by only aiming at reduced information. Namely, 
we consider a partition A = {Ai, A 2 ,..., A^} of S (with m ^ N). Now sample m individuals 
from the present population and follow back the ancestry of the sites in Ai in the first individual, 
in A 2 in the second individual, ..., in Am in the m’th individual, without considering any other 
sites and any other individuals, as in Figure [3l That is, each locus is considered in one individual 
only. The result may be viewed as a marginalised version of the ancestral recombination process, 
and, in the diffusion limit, turns into a marginal version of the multilocus ARG starting from a 
sample of size m. We will see that this information is sufficient to characterise the time evolution 
of the expected linkage disequilibria of all orders. We will not employ any scaling or limit, in 
order to allow for arbi trary strengths of recombination. It will turn out that the approach of 
( 20101 ) actually corresponds to this marginal ancestral recombination process, 


Bobrowski et al 


although this is not apparent from their formulation forward in time. 

More precisely, the letters at the loci considered at present, together with their ancestry, can 
be constructed by a three-step procedure (see Figure S]). First we run a partitioning process 
on IP(5), backward in time, starting at a given initial partition Uq with iJCol = m. Et 
describes the partitioning of sites into parental individuals at time f; sites in the same block 
(in different blocks) belong to the same (to different) individuals. Glearly, \Et\ is the number 
of ancestral individuals at time t. The process {Et}t^o is independent of the types and will 
be described in detail in the next section. In the second step, a letter is assigned to each site 
of S at time t (i.e. in the past) in the following way. For every part of Et, pick an individual 
from the initial population (without replacement) and copy its letters to the sites in the block 
considered. For illustration, also assign a colour to each block, thus indicating different parental 
individuals. In the last step, the letters and colours are propagated downwards (i.e. forward in 
time) according to the realisation of i Et}t>n laid down in the first s tep. A similar construction 
was used in the ancestral process by IBaake and von WangenheimI (j2014l ). but restricted to a 
sample of size 1 (i.e. start with lio = 1), and in discrete time in the deterministic limit. Let us 
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now describe the partitioning process in detail. 




Figure 4: Construction of one possible ancestry of a collection of sites that correspond to the initial partition 
Ao = {{!}, {2, 4}, {3, 5}}. The upper panel shows the partitioning process (backwards in time). In the lower 
panel, letters and colours are assigned to each block of Et and propagated downwards (forward in time). 


4 The partitioning process 

The partitioning process {X!t}t^o is a Markov process on ^(S), which describes how the sites are 
partitioned into different individuals backward in time. Since there is a one-to-one relationship 
between the individuals and the blocks of the partition, we may identify individuals with the 
ancestral material they carry. 

The process consists of a mixture of splitting (S) and coalescence (C) events. It can 

be constructed independently of the types. In this section, we describe the process by arguing 
on the grounds of the underlying Moran model; in Section [Gj we will formally prove that this is 
indeed the correct dual process for it. 

Since we trace back sites in subsets U E S (rather than complete sequences), we need the 
corresponding marginal recombination probabilities 

■= Yj ^A ( 9 ) 

.4eo^2(S') 

A\jj=B 

for any B G 0<g2(t^), where r^ = r^. Note that, for \U\ = 1, the only recombination parameter 
is r^ = 1. If [/ is ordered in S (i.e. U = {x € S : min(C/) ^ x ^ max([/)}) and B / l|fj, then 
r^ is simply the probability of crossover after the (unique) site that leads to partition B. If U is 
not ordered in S, then r^ is the sum of the probabilities of all crossovers that lead to partition 
B, as illustrated in Figure [5l 

Assume now that U is an unordered block of Ut- This means there is so-called trapped material, 
that is, non-ancestral sites enclosed between ancestral regions. All crossover events within a 
given trapped segment contribute to the separation of the adjacent ancestral segments - in 
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contrast to crossovers in flanking non-ancestral regions to the left or the right of U, which do 
not affect the genealogy. Note hnally that the upper index in can, in principle, be omitted 
since U = Bi, and we will do so when appropriate. 




+ Q- 


+ Q- 


Figure 5: Let S = 5} and U = {1,4,5} C S. For the partition B = {{1},{4, 5}}, there are three 

recombination events that partition U into B, thus { 2 , 3 , 4 , 5 }} + r{{i, 2 },{ 3 , 4 , 5 }} + r{{i, 2 . 3 },{ 4 , 5 }}. 


Now start with the initial partition Eq. Suppose that the current state is Et = A = ..., Am} 

and denote by A the waiting time to the next event. A is exponentially distributed with 
parameter m, since each block corresponds to an individual, and each individual is independently 
affected at rate 1. When the bell rings, choose a block uniformly. If Aj is picked, then A’t+A is 
obtained as follows (see Figure E] for an example). 

In the splitting step, block Aj turns into an intermediate state A'- = J with probability 
J € 0^2 (^j)- In detail: 

(51) With probability the block Aj remains unchanged. The resulting intermediate 
state (of this block) is A'- = Aj. 

(5 2 ) With probability J € 02(Aj), block Aj splits into two parts. A'- = J = {Aj^, Aj^}, 

which are ordered in Aj, but not necessarily in S. Recall that, via ([9]), takes into 
account all recombination probabilities that lead to J, including those within trapped 
material. 


Now, each block of J chooses out of N parents, uniformly and with replacement. Among 
these, there are m — 1 parents that carry one block of AM\j = A \ Aj each; the remaining 
N — [m — 1) parents are empty, that is, they do not carry ancestral material available for 
coalescence. Coalescence happens if the choosing block picks a parent that carries ancestral 
material; otherwise, the choosing block becomes an ancestral block of its own, which is available 
for coalescence from then onwards. The possible outcomes are certain coarsenings of AM\j U J, 
namely: 


J = {Aj} (case (Si)), then either 


(Cip) With probability [N — (m— l))/iV, block Aj does not coalesce with any block of AM\j- 
As a result, Et+A = Et = A. 

(Cl, 2 ) With probability 1/A, block Aj coalesces with block A^, k & M \ j. This results in 
Et+A = -4M\{j,fc} 


If A = {Ajj, Ajj} (case (S 2 )), we get the following possibilities: 


(C 2 ,i) With probability {N — (m — 1))(A — m)/N‘^, no block of J coalesces with a block of 
AM\j , so Et+A = AM\j U J. 

(C2,2) With probability {N — (m — 1 ))/A^, one block of J coalesces with block A^, k E M\j, 
while the other block of A chooses an empty individual. This ends up in the state 
Et+A C Aj 2 } or Et+A C {Ajj^ ^jj-, Ajj}. That is, in going 

from Et to Et+Ai either block Aj^ or Aj.^ is moved from Aj to A^. 
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(C 2 , 3 ) With probability {N — (m — 1))/N‘^, the blocks Aj^ and Aj^ coalesce with each other, 
but choose an empty individual, which gives Ut+A = A. 

( 02 , 4 ) With probability the block Aj-^ coalesces with A^ and Aj^ coalesces with A^, 

k,£ E M\j. This yields either St+A = AM\{j,k,e} U if fc / £, or 

^t+A = -4M\{j,A:} if ^ 


o- 

IZ3 




1 (^ 2 , 2 ) 



^t+A 


{St \A2)yjj 

St 


Figure 6: One step of the partitioning process with cuexitrrent state St = {Ai, A 2 , A 3 } = {{!}, {2,4}, {3}}. In 
this example, A 2 is chosen and splits into J' — {{2}, {4}}. In the following step ( 02 , 2 )! the leading part coalesces 
with Ai, whereas the trailing part remains separate, so that we end up in St+A ~ {{1, 2}, {3}, {4}}. 


Summarising, we see that a transition from A to B, via partitioning of block Aj into j E M, 
J E 0 !^2{Aj), is possible whenever B AM\j U J' and = AM\j, of, equivalently, 

whenever 

B\a.>J and B\a^^. = AM\j- 

Each block of coalesces into every block currently available with probability 1/N, and remains 
separate with probability {N — k)/N if there are currently k blocks available; in the latter case, 
the block considered becomes number fc + 1. We can therefore summarise the rate of the said 
transition as 




( Aj 1 (jv-(m-l)) 

I W TVI^T (Af-|B|)! 

lo, 


if B\a^ > J, B\a^^^ 
otherwise. 


A 


M\j, 


( 10 ) 


Note that this includes silent events where B = A. Thus, the partitioning process {St}t'^o is a 
continuous-time Markov chain on P(S') characterised by the generator Q := (Oab)A,B&¥{S) with 
nondiagonal elements 


^AB = X] X] 

j(zM J'G0^2{Aj) 

Aj 1 (jV-(m-l))! 
'j W (Af-|B|)! ’ 


= 


2 

W 




S _|_ At 




if ^\Aj = = AM\j, for some j e M, J £ 02{Aj), 

if B = AM\{j,k} U for some j ^ k £ M, 


0 , 


for all other B A A, 


( 11 ) 


and = ~ J2bgP(S)\A^AB- Note that, for J" E 02 (Aj) we have distinguished between 

B\j<^. = J and B\j<^. = 1|^. J. The latter corresponds to fc = f in (€ 2 , 4 ) and leads to the 
same transition as a pure coalescence event in (Ci^ 2 ) More precisely, the total coalescence rate 
of j and k is 


1 

N 


(rf^ +ri'“) + 


1 

iv2 V 


( 


E 

j€ 02 iAj) 


r^-' + r 


Ak 

K 


IC£ 02 {Ak) 


2 

w 






as stated, since X)j'G 02 (r/) ~ ^ ~ U S S. Note also that transitions to partitions B with 

\B\ > N are impossible, as it must be. 
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Rema rk 3. In fact, this gene r ator 0 coincides with the generator 0 worked out bv lBobrowski and Kimmel 
( 20031 ) and Bobrowski et ah ( 2nifll i with a very different approach, forward in time. For n ^ 

3, they state the generator matrices explicitly, and the identity with m is easily checked 
by elementwise comparison. For n > 3, they provide an algorithm, which runs through all 
individuals and all sites and builds up the matrix 0 incrementally, in the following manner. For 
every given individual, leading and trailing segments (for the split after site i, for all i G S\n) are 
taken into account, irrespective of whether the segments contain ancestral material. This way, 
the algorithm does not distinguish between transitions induced by recombination events within 
ancestral (or trapped) material and recombination events that are invisible in the genealogical 
perspective, that is, events that are effectively pure coalescence events. Instead, a case distinction 
is performed that is based on whether or not one or both segments coalesce with individuals 
that do or do not carry ancestral material. A detailed investigation of this approach, which 
involves expanding the cases into 11 subcases and rearranging these according to the emerging 
partitions of the complete ancestral material, leads precisely to our cases (C 2 ,i) to ( 02 , 4 ) (here, 


both emerging segments contain ancestral material) and (Ci^i) and (Ci^ 2 ) (here one segment 


is empty). Since this approach somehow disguises or mixes the various partitions of ancestral 
material that may arise due to a transition, it does not lead to a closed expression for 0. In 
contrast, our approach yields the matrix elements explicitly for arbitrary n, and gives them a 
natural and plausible meaning in terms of the partitioning process in backward time. <0> 


Limits of the partitioning process. We now examine how the partitioning process behaves in 
the two limiting cases mentioned in Section [3l namely, the deterministic limit and the diffusion 
limit. Recall that, in the deterministic limit, we let —>■ 00 without rescaling the recombination 
probabilities or time. Consider, therefore, the family of processes N = 1,2,..., 

generated by 0 ^^\ where we again make the dependence on population size explicit through 
the upper index. In the limit, only the pure splitting events (C 2 ,i) survive, more precisely: 


Proposition 1 (Deterministic limit). In the deterministic limit, the sequence of partitioning 
processes with initial states = cr converges in distribution to the process 

with initial state Uq = a and generator 0' defined by its nondiagonal elements 


^AB — 



if B = -^M\j ^ some j & M and J € 02(Aj), 

for all other B A. 


Hence, is a process of progressive refinements, that is, ^ for all t > t. In 

particular, if Hq € 0{S), then G 0(5) for all times. 


Proof. Inspecting the ^dependence of the elements of 0 = 0^^'l in Eq. m gives the following 
order of magnitude for the nondiagonal elements: 


0im _ 


Amw '■j (1 + O(ir)), if = J. tor j e M , J e 02(.4,), 

+ r’t'') + if ^ ^ if fffj.ii} fiii' i # fi e V, 

0, for all other B / A. 


( 12 ) 

Obviously, 0^^'l = 0' + 0{1/N), which proves convergence of the sequence of generators 
of to that of This entails conver g ence of the corresponding sequence of 

semigroups. By Theorem 4.9.10 of lEthier and Kurtzi (jl985l h this guarantees convergence of 
to in distribution. 
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The remainder of the statement is obvious since, under 0', the only transitions are those that 
involve the refinement of a single block, say Aj, into two blocks ordered in Aj. If Sq is ordered in 
5, then all its blocks are ordered in S', and all blocks of Et will be ordered in S for all times. □ 


Remark 4. Obviously, in this limit, there are no coalescence events, so that ancestral material 
that has been separated once will never come together again in one individual. When starting 
with Eq = {S}, the genealogy may be represented by a binary tree, which successively branches 
into smaller segments; for other initial conditions, one gets a corresponding collection (i.e. a 
forest) of binary trees. We call these trees ancestral rec o mbin ation trees or ARTs\ a discrete-time 
analogue was studied by Baake and von Wangenheim ( 20141 1. ‘C> 


We now turn to the diffusion limit and use the factor N rather than (the more common) 2N 
since our N is the haploid population size. Here, one considers a sequence of processes in which 
time is sped up by a factor of N and the recombination probabilities are rescaled such that 
limAT^oo Qai 6a ^ constant, for A G ©2(5'); consequently, —>■ 1 as —>■ oo. Note that 

the are rates rather than probabilities. The corresponding ARG is the obvious generalisation 
of Hudson’s original ARG to n loci, which we formulate here in our framework for the sake of 
completeness, as follows. Every ordered pair of lines coalesces at rate 1; every line splits into 
two at rate for every A G ©2(5'), and the ancestral material is distributed between the new 
lines according to A. 


In this formulation, however, certain silent events are included, namely those events that happen 
in non-ancestral material flanking the ancestral parts. These events do not affect the partitioning 
of ancestral material and may be removed by working with the marginalised recombination rates 
instead. That is, if a sequence currently carries a set [/ of ancestral sites, then the relevant 
recombination rates (in the diffusion limit) are with B G © 2 ( 1 /), which are defined as in ([U]) 
but with r replaced by g. If we now restrict attention to the ancestry of n loci partitioned 
between m individuals, we obtain the marginal version of the ARG, which may be formulated 
as follows. 


Definition 2 (Marginalised n-locus ARG). Start with the set of n sites distributed across m ^ n 
individuals (or lines) according to a partition Eq with m parts. Throughout the process, every 
line is identified with the ancestral material it carries. If it currently carries ancestral sites 
[/ C S', it splits into G ©2 at rate gj. Every ordered pair of lines coalesces at rate 1, and 
so do the ancestral sites they carry. That is, the marginalised ARG is the partitioning process 
defined by the generator 0 " with nondiagonal elements 


0.4B = <: 


Aj 

ej 

2 , 

0 , 


if H = AM\j U J for some j G M,J G © 2 (Aj), 
AB)- A and \B\ = |.4| - 1, 
for all other B A 


Proposition 2 (Diffusion limit of the partitioning process). In the diffusion limit, the sequence 
of partitioning proeesses with initial states E^^ = a converges in distribution to the 

process {E'f}t^Q with initial state E'f = a and generator 0". 


Proof. Due to the rescaling of time, the generator of {E^J}t^Q has nondiagonal elements N0^ ■ 
Referring back to (fT^ . they converge to limAr_).oo = 0_4^, since we have ^ 1 and 

Nrj -A gj for ff G 02(17)■ With the same argument as in the proof of Proposition [H one 
obtains convergence in distribution as claimed. □ 


As was to be expected, only pure splitting events and pure coalescence events survive in the 
diffusion limit. The ‘mixed transitions’, which involve both splitting and coa lescence (i.e. the 
dashed lines in Figure [2]) vanish under the rescaling; see also Hein et al. ( 20051 . Fig. 5.11). 
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5 Recombinators and sampling functions 


In this section, we will have a closer look at three operators associated with recombination and 
how they are related to each other. We start by generalising our recombinators, then introduce 
closely related sampling functions and finally multilocus correlation functions, known as linkage 
disequilibria. 


Rocombinators. We have already met for A G 0^2(5'); we now need the generalisation 
to arbitrary A G P(5). For u G AI+(X) \0, we first define the non-normalised recombinator via 

(8) • • • 0 :, (13) 


where the symbol : ... : means ‘site ordering’ and is used to indicate that the product measure 
refers to the ordering of the sites as specified by the set S. When the meaning is clear, we drop 
these extra symbols. In words, turns a; into the product of its marginals with respect to 
the blocks in A. We will throughout denote non-normalised mappings by an overbar. Clearly, 
R 0 {uj) = ||a;||, Ri{u}) = ui and ||ii_ 4 (w)|| = The corresponding normalised version is 


-Ryt(w) := 


Ra{uj) 


(14) 


which is well-defined since w 7 ^ 0. Obviously, Rji^{u) = i? 4 (a;/||a;||) and Rji^{oj) is a probability 
measure on X, which coincides with ([ 6 |) for A G 0^2 (*S'). 

Let us now give a probabilistic interpretation for the case that a recombinator R^ acts on a 
certain population described by a counting measure z G E. For the moment, attach labels 
l,2,...,iV to the N individuals in the population, and let these individuals have (random) 
types ... ,X^ G X at time t. The type distribution then is Zt = For U C S 

and k G {1, • • • ,X}, let X^^ := 'Kjj{X^), and consider the following procedure. Let a partition 
A = {Ai,..., Am} of S together with a collection of labels i = {£ 1 ,. ■ ■ Am) G {1,..., A/'}'” 
associated with the blocks be given, i.e., (A,£) ■= ■ ■ ■ 1 {-^mAm)}- Then, piece together 

a sequence by taking the sites in Ai from individual £ 1 , the sites in A 2 from individual £ 2 ^ ■ ■ ■ 
the sites in Am from individual £m- The resulting sequence is Xfj^ := We 

are now interested in the event 


{ W ,^ = x }:= U {xIj, = x] 

and the corresponding counting measure = 3;}|. Clearly, this counts how often one obtains 

sequence x when performing the above procedure on a population Zt and all combinations of 
individuals are included. Let us emphasise that individuals are combined with replacement, that 
is, two or more blocks may come from the same individual. Therefore, the event {W,.A = x} 
may also be understood as the union of the independent events {Xt A' = Xa.}, j G M, where 


U ( 15 ) 

Therefore, 

\{Xt,A = a:}| = n =XAi]\= n {'^A^-Zt){xA^) = {RAi^t)){x). (16) 

j&M j£M 
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Clearly, Rj^[Zt), the corresponding normalised version, is the type distribution that results when 
a sequence is created by taking the letters for the blocks in A from individuals drawn uniformly 
and with replacement from the population Zf. So 

{Ra^^)) (®) = P [^t,A = x\Zt = z], 

where P denotes probability. Note that the left-hand side depends on time only through the 
value of Zt- 

Sampling function. For A € P(5) and OJ € A^ + (X) \ 0, we now define our sampling function 

(17) 

B)pA 

where is the Mobius function in Eq. ([2]). is not a positive measure in general; but it 

will turn out as positive for the important case where w € E with ||a;|| ^ \A\, see Lemma [T] 
We will therefore postpone the normalisation step. In any case, Mobius inversion (compare ([3]) 
and Q) immediately yields the inverse of (HZl): 

Fact 1. For every A G P(5'); 

Rj,{u:) = Y, □ 

B)pA 

We can now give a meaning by reconsidering the procedure that led to (1161) but, this time, 
individuals are not replaced. That is, for |M| ^ N, we now look at the events 

{W,^ = x}:= U {xIa = ^] (18) 

AT}"* 

and the corresponding counting measure |{W,yl = 3:}|. Since individuals are not replaced, the 
events {Xt a = x a.}, j G M (defined as in (|15p with X replaced by X) are now dependent', an 

expression for \{Xt^A = x}\ analogous to (1161) is therefore not immediate. Instead, we resort to 
an inclusion-exclusion argument and prove 

Proposition 3. For A G P(5) with |.4.| ^ N and Zt G E, we have 

\{YA = x}\ = {H_^{Zt)){x). 

Proof. Fix a given partition A G P(5) with |^| = m ^ N. For every i G {1,2,... , iV}”*, the 
pair (A,£) uniquely defines a pair where £ G G {1, 2,... , ^ ik'^ j A k} 

and ^ as follows. Join all blocks of A that have the same label, and attach that label 
to the new block. The result is The other way round, every {B,i) with B A and 

G G {1, 2,... : ij A j ^ k} uniquely dehnes the labelling i of the blocks of A 

(keep in mind that A is fixed): block Ak € A receives the label of that block Bj G in which 
it is contained. We can therefore identify the set {{A,i) : i G {1, 2,... , A^}”^} with the set 
[j , i) : i G ^ {1,2 ,... : £j A j 7^ k}. With (fTHI) and (fT8|) in mind, we can 

therefore decompose the event {W,.A = x} = which entails 

B^A 

By (jl6l) . the left-hand side equals [R_^{Zt)){x). Due to the Mobius inversion principle (applied 
backwards), \{Xt^B = 3:}| on the right-hand side must equal (F[j^(Zt)){x), as claimed. □ 
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Lemma 1. For A G P(5') with |^| = m ^ N and z a E, H_^{z) is a positive measure with 

||F 4 (z)|| = N {N - 1) ■ ■ ■ {N - m + 1) > 0. 

Proof. Since, under the given assumptions, {H_^{z)){x) = \{Xt^_A = x \ Zt = z}\ ^ Q for all x by 
Proposition [3l it is a positive measure, and its norm can be evaluated via 

\\M4\=Y.\{^hA = x\Zt = z]\. 

By means of (I18p , this equals the number of possibilities of how to choose m labelled individuals 
out of N ones without replacement, where the order is respected; this is N (A^ —1) • • • (N — m+1), 
which is positive since m ^ N. □ 


Under the assumptions of Proposition [3l we can therefore define the normalised version of H_^{z)'. 


H^{z) := ^ 




Ha^z) 




(19) 


Ha{z) is the type distribution that results when a sequence is created by taking the letters for 
the blocks as encoded by A from individuals drawn uniformly and without replacement from the 
population z, hence 

{HA{z)){x)=F[Xt,A = ^\Zt = z]. 

Hj^ will later serve as duality function. The situation described here is exactly what happens 
when a sample is taken in our marginal ancestral recombination process: either the initial sample 
(according to Uq, from the present population Zt) or the ancestral one (according to Xt, from 
the initial population Zq) - hence our name sampling function. In this light. Fact [T] expresses 
counting with replacement in terms of counting without replacement, provided a; is a counting 
measure. 

It is also instructive to express the normalised sampling functions in terms of the normalised 
recombinators. For z & E and |^| ^ N, this gives, via (jl4l) . 

B)pA 


Note that {N — |.4|)! This illustrates how the inclusion of coarser 

partitions yields higher-order correction terms. The other way round, using (mD, Fact m 
and (fTUl) . one gets 


HAV = E 


A^! 




An^i(Ar- |B|)! 


H^{z). 


( 20 ) 


Restriction to subsystems. Recall that we write the restriction of a measure tv G Ad+(X) 
to a subspace X[/ := of X as tTjj.oj := lo o which corresponds to marginalisation. 

Clearly we can also define recombinators for any non-empty subset U O S and any partition 
A = {^ 1 ,..., Am} G F{U) as Ra{ttij.uj), in perfect analogy with Ra{^) for A G F(S), which 
is Ra{w)] and likewise for Ra, , and Ha (if a; 7 ^ 0). For clarity, we sometimes denote the 
subsystem by a superscript. However, as in the case of the marginal recombination probabilities, 
it can be dispensed with since U = Aj if ^ G P(C/). The interpretation in terms of sampling, 
as well as Fact [H carry over. 

Let us collect some basic properties of recombinators: 
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Fact 2. For A,B £ P(S') and U,V C S with S = U OV one has 


fAJ — R^ab- 


(B) TTjj.R^iu) = R^^^iTTu-Uj). 

(C) If in addition A ^ {[/, F}, then R^ = Ra\^ ® ^^\v' this reads 

R^{uj) = 0 RA\y){^) = {Ra\u^'^u-^)) ® {RA\yiAv^))- 


Here and in what follows, we may omit the argnment when the meaning is clear. 


Proof of Fact^ Property |(A)| is Proposition 2 and property (B) is Lemma 1 of iBaake et al 


( 2 OI 5 I I (they both remain true in our normalisati on). Property (C) is an obvious generalisation 
of Proposition 2 of Ivon Wangenheim et ahl ( 20101 1. It is easily seen by using first property |(A)[ 
then (1131) , then |(B)| and finally (1131) once more to give 


Rjico) = R^u,v}{RXi^)) = {{^u-RX) ® {7Ty.RX)){u;) 


= {RA\ui'^U-‘^)) ® = {Ra\^ ® RA\y)i‘^) 




■F>U 




□ 


Let us note a connection between recombination and sampling that will be important in what 
follows. 

Lemma 2. Let S = U (JV for two nonempty subsets U,V C S'. For two partitions A G P{U), 
B G P(P), the recomhinator and the sampling operator satisfy 

C^AUB 

C\y=B 

Proof. Using (fT71) followed by Fact [2] [(^ and Fact [1] we get 

rX®h^ = rX(^[Y .= E= EE 

^ V)aB V)pB V^pB E)pVliA 


Changing the summation order and applying ([2|) finally leads to 

Ra®Hb= E E E ^c- □ 

OpAvJB B^V^C\y C-pA^B 

C\y=B 

Remark 5. In a perfectly analogous way, one can show 

hUHb= E Hi 

OpAvJB 

C\^=AC\y=B 

This illustrates once more that, unlike the R 4 , the do not have a product structure; this 
reflects the dependence inherent to drawing without replacement. < 0 > 
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Correlations (or linkage disequilibria). Linkage disequilibria (LDE) are used in population 
genetics to quantify the deviation from independence of allele frequencies at the various sites in a 
sequence. From thre e sites onward s, many different notions of linkage disequilibria are available 
in the literature, see iBiirgerl ([20001, Chap. V.4.2) for an overview. 


We will use as LDEs the ge neral c o rrelat ion functions, which are widely used in statistical 
physics, see Dvson ( 1962l l or Mehta (1991, Chap. 5.1.1). This results in an explicit formula 


for multilocus LDEs for an arbitrary numbe r of sites in ter ms of sums of products of mar ginal 
frequencies, see also Baake and Baake ( 2003 . Appendix) or Gorelick and Laubichler ( 2004l l . As 
we will see, common definitions for two and three sites coincide with ours. 

For any given subset [/ C 5 and A € F([/), we first define correlation operators as 




( 21 ) 


B^A 


Note that the summation is now over all refinements of A, in contrast to our sampling functions, 
which involve all coarsenings of A. The restriction to subsystems stems from the fact that one 
usually considers deviation from independence on small subsets of S. 

The have a product structure, > which is obvious from (|2ip together with the 

product structure of the recombinators (Fact E] [(C)] ) and that of the M5bius function (Eq. ([ 2 ])). 
Eq. (|21l) has the inverse 

\B\ 

B^A B^A j=l 

due to inversion from below (see Section [2]) . The latter can be reformulated as 

\B\ 

La = Ra-Y. ( 22 ) 

B^A j=l 


The case A = 1|;7, U C S, now is of special interest. In line with population-genetics 
understanding, we define the multiloeus linkage disequilibrium with respect to the sites in U 
by letting act on the marginal measure ttjj.uj, uj G M_|_(X) \ 0: 

L^{7Tu.w)= p,{A,l\u) ^Ai'^U-^)^ 

A&iU) 

cf. Eq. (|2ip . Note that (ttjj.lv) is again a measure on 7 ru(X), but no longer positive in general. 
With the help of (1221) . it can be reformulated as 

|B| 

L^inu-uj) = R^iTTjj.uj) - X] n 

B^i\u i=i 


which is Eq. (1) in iGorelick and Laubichlen (12(10411. Likewise, t his alternative form ulation of 
multil ocus LDEs agrees with previous ones from Geiringer ( 1944l l . Bennett ( 1954l l. and Hastings 
(1984) up to \U\ ^ 3. 

Example 2. For S = {1,2,3,4} the LDE with respect to the sites in U = {2,4} reads 



.w)(x) 


/, n /, M 
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0 

The correlation operators can also be expressed in terms of our sampling operators. Eqns. m 
and ()2np . together with a change of the summation order, lead to 


- 




m 


I3=4A 




(A^- |C|)!iViB| 


Hr^ = 


2^ He 

C&{U) 


u 


E 

B4A/\C 


m 


(A^- |C|)!iViB| 


KB. A). (23) 


For a counting measure z € E and U C S with |f7| = A; ^ 3 ^ A^, Eq. (1231) yields a particularly 
nice explicit expression for the LDEs: 

LiiT^uK = j.^k(N'-kV ^ KA,1\u)HaKuK, (24) 

^ A£P{U) 


as is easily verified. For larger k, the explicit formula gets more involved. 

Let us now consider for A € P(t/) \l|f/. Due to its product structure, the collection of all 
Li {TTy.u), V C U, determines all A € P(C/). This is why, for a deterministic ui, the 

Lj^Ku-uj). A A 4|[/) of iio particular interest of their own. This changes, however, when oj is 
random (like Zt). For we typically do not know the law of Zt completely; rather, we have access 
to the expectation of certain functions of Zt- More precisely, let (/? be the law of Zt and 
denote the expectation with respect to ip (that is, for a function / of Zt. E^[/] = f f{z) (iip{z)). 
It is important to note that the product structure of the recombined measure does not carry over 
to the expectation. That is, for A G P(D), -Zt)] A Ra^^^v Ajj.Zt]) in general, see the 

discussion in iBaake and HermsI (j^OOslL Conseouentlv. E,.,[L^('7r^r.Z/')] A ^ .E,., [Zt]) 

in general. In the stochastic case, therefore, it is interesting to consider the for A A 


as well. The expectations EKl^Ku-Zt)] contain information on how the mean LDEs 
part of the sequence depend on the mean LDEs in other parts of the sequence. In the next 

section, we will obtain an ODE system for the E(^[L/^ {Zt)]. A G P(5'), and these translate into 
s s 

E(p[i? 4 (Zt)] and thus into E(^[L(^(Zt)] via (12311 . Marginalisation can then be used to calculate 
the corresponding quantities on [/ C 5, such as E<^[L^(7rfj.Zt)] for A G P(C/). 


m one 


6 Duality 


Duality is a powerful tool to obtain information about one process by studying another, the dual 
process. The latter may, in an optimal case, have a much smaller state space than the original 
one. Duality results are essential in interacting particle systems in physics and in population 
genetics. They are often related to time reversal. The most famous example in population 
genetics is arguably the moment duality between the Wright-Fisher diff usion forward i n time 
and t he block c ounti ng process of Kingman’s coalescent backward in time ( Donnelly 19861 : Mohle 
2001). Mano ( 2013l l extended this result by incorporating recombination into the two-locus, 
two-allele case. His results are based on the original version of the ARG and thus on the 
diffusion limit. 

We will briefly explain the general duality concept and then prove that our processes {Zt\t^Q 
and {Et}t^o are duals of each other. For the general principle, let X = {Xt}t^o and Y = 
be two Markov processes with state spaces E and F. Define hy M (Ex F)t, the set of all bounded 
measurable f unctions on E x F. The following definition of duality with respect to a function 
goes back to Ligefett ( 19851 ): see also the recent review by Jansen and Kurt (2014). 
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Definition 3 (Duality). The Markov processes X and V, with laws (p and ip, respectively, are 
said to be dual with respect to a function H M{E x F)b if, for all x £ E, y ^ F and t ^ 0, 

[H{Xt,y) \Xo = x]=E^ [H{x, Yt) \Yo = y]. (25) 


If E and F are finite, every function H £ M{E x F)b may be represented by a matrix with 
bounded entries H{v,w), v £ E, w £ F. If, further, X and Y are time-homogeneous with 
generator matrices 0 and A respectively, the expectations in equation (|25p may be written in 
terms of the corresponding semigroups, i.e., 

[H{Xt,y) I Xo = x] = ^(e*®),,if(u,y), 

(26) 

E^ [H{x,Yt) \Yo = y] = (e*^)^^il'(x, tc). 

wGF 


Since the duality equation (j25p is automatically satisfied at t = 0, it is sufficient to check the 
identity of the derivatives at t = 0. That is, Eq. (j25D holds for all times if and only if 


-E^ [HiXb,y)\Xo = x]l^^ 


'^A^yH{v,y) 

vGE 


H{x,w) 0yw = 

w£F 


dt 


E^ [H{x,Yt)\Yo = y]l^^ 


(27) 


for all X £ E, y £ F. As a short-hand of Eq. (1271) . one can write AH = H0^, where T denotes 
transpose. 

We will now present a duality result that justifies our construction of a marginalised sample 
at present via the partitioning process and sampling from the initial population (cf. Figure 0]). 
Indeed, it is not coincidence that we have denoted our sampling functions by and our 
generators by A and 0. 

Theorem 1. The population process {Zt}t^o the partitioning process {Xt}t^o with the 
generators A and 0 and resulting laws p and respectively, are dual with respect to the sampling 
function H defined in Eq. (HSl). Explicitly, 


[H_A{Zt) I Zo = z] = E^ I Ao = M] (28) 

for all A £ P(S) and z £ E. 


Before we embark on the proof, let us briefly comment on the meaning of this result. 

Remark 6. Eq. (j28p is the formal equivalent of the construction in Figure 01 To see this, 
recall the random variables from (1181) . With their help, the left-hand side of ()28p may be 
reformulated as a probability distribution, 

\_H_A{Zt) I Zq = z] = E(^[P = •] I Zt, Zq = z] = Py, = • I Zq = zj, 


since the expectation is o ver all realisation s of Z ^. 
distribution considered by Bobrowski et al. ( 2r)l(ll l. 
equal to 


The right-hand side is the probability 
Likewise, the right-hand side of (1281) is 


E,/, [Hj]fiz) I Xo = ^] =E,p [P[Xo,i7t = •] I Xq = A\ = P.^ [Xo^Et = • I = -4] , 


since the expectation is over all realisations of Ft. The right-hand side is the distribution of types 
when sampling from the initial population according to the partition Ft, where it is understood 
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that the initial population consists of the types Xq, ... ,X^ with S^k = z. Recall that 

time runs forward in Zt, X^, and Xt^_Aj but backwards in X*. <C> 

In order to avoid case distinctions in the calculations in the remainder of this section, let us agree 
on the following conventions concerning (partitions of) empty sets. Namely, we set A 0 := 0, 
= R 0 {tt 0 .z) := TT^.z = Ijzll = N, and /i(0, 0) := 1. We now collect some auxiliary 
results in the following Lemma. 

Lemma 3. Consider a counting measure z & E, a partition A G P(S') with |.A| = m ^ N 
and corresponding index set M = {1,... ,m}, and a partition B € P(S'). Then, the following 
statements hold: 

(A) Y,{Rb{z)){x) [H^{z + 5^)-H^{z)] = 

xex j&M ^ 

(R) - ^y) - = -mH_^{z). 

ygx 


Note that the left-hand side of statement (B) is always well-dehned, since z — 6y < 0 can only 
occur with z{y) = 0, in which case the term vanishes. Note also that, with the above conventions, 
the right-hand side of identity |(A)| evaluates to 


j&M " 


Proof. We first observe that 

^ {RBiz)){x){Tru.6,,) = TTJJ.J2 {RBi^))ix) s,, = Tru.{Ri^{z)) = R^i^iTr^.z) (29) 

icEX icEX 

by Fact [2j We next evaluate X^xex (-^b(-2))(^) [^Ai^ + — Ra^^)] by expanding R 4 to 

separate the action on z from that on 5x, summing against Rjs{z) (using (| 2 ^ L applying Fact[T] 
and changing summation: 

E (^b(^)) {x) [Rj^{z + 5x) - Ra{z)\ 

31GX 

xex 0^jcM 

= E E E {Bc-sr.iJa) 

0^JCM 0^J(ZM C)pA^,^j 

\V\ 

= E E ® Rb\^^{^)^ 

V)pA j=l ^ 


where, in the last step, every Aj reappears as one Dj. Using this together with (I17D and ([5]), 
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we obtain 


Y1 {^Biz)){x) [H^iz + 5cc) - h{A,C)Y {RBiz))ix) [Rdz + 5^) - Rdz)] 

xex C^A xeX 

/- X 

= 'Y X/ X/ ® ^B\ij .) (^) 

C-^A V-^C j=l ^ 

\V\ 

= '^ '^{YT^\Dj® Y ='^ {^AM\j® ^B\^'){z)^ 

V)pA 3=1 ^ A^C^V j&M ^ 


which is statement |(A)[ 

In an analoguous way, we can prove statement |(B)[ 

Y^^y^ - ^y) - Ra{^)\ = Y (-1)''^' {RAmY^Am\j-^)) ® {Rt'iT^Aj-^y)) 

y€X 0^JCM y€X 

0jlJCM 0^JCM 


= E (-1)'''' 

0jlJCM 


|C| 


E = E ^c(A E 

B'^Am\j'JAj C^A j=l 


E 

0jlKCCj 


\C\ 

= E^c(^)E [(i-i)'^^'-i 


= -Y\C\Hciz), 


C^A j=l C)pA 

where, in the second-last step, every Aj reappears as a Cj. We therefore get 

Y - 6y) - H_^{z)] = Y Y1 i^Biz - 6y) - 

yex B>pA j/gx 



= -Yf^i-^^>3)Y\C\Hciz) = -Y\^\^ciz) Y 

B^A CVB CvA A^B^C 

= -\A\H_^{z), 


as claimed. 


□ 


We can now proceed as follows. 


Proof of Theorem^ We start with the partitioning process. We first note that 


E 




(^N — [m — 1))! 

(A^-I^l)! 


N\x\ 


(30) 


for j ^ M and \J'\ ^ 2. This is easily verified by direct calculation; namely, for \J'\ = 1, the 
sum on the left-hand side equals [N — {m — 1)) -|- (m — 1) = N; for \ J'\ = 2, it evaluates to 

[N — {m — 1)) {N — m) + {N — (m — 1)) (2m — 1) -|- (m — 1)^ = N"^. 
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We now use the formulation of the process via (jlOp and (jlll) in the first step, normalisation 
and (|30p in the second, Lemma [2] in the third, and finally another normalisation step to calculate 


^AB 

BeP(S) 


= 


E E 


'J 


-{m- 1))! , 


ivi-71 ^ 




E E 

j&M J'&0^2(Aj) 


'J 

N\J\ 


E I - ffA I (•-) 




/GAf j7'GO<2(24i) \ / 


j£M J G0^2(^j) 

E E 

jGM J'£0^2iAj) 


(31) 


We now turn to the type distribution process. Here we first evaluate, with Lemma [3] (B), that 

E + 4 - Sy) - H_^{z)] 

= '^{z + 4)(y) H^iiz + 4) - 6y) - '^{z + 4)(y) Hj^{z) 


— E^^ ^a;)(y) + 5x) — 6y) — Hj^{z + 6x) + Hj^{z + 5x) — Hj^{z)\ 

yex 

= {N + l-m) [H_^{z + Sx) - H_^{z)] - m H^{z). 

From this, we obtain via summation against Rq{z) and use of Lemma [3| | (A) | that 

E E (®) ^(y) [^Ai^ + 4 - Sy) - H_^{z)] 


xSXj/GX 


= (iV + l-m) ^ ^ ^B\^)iz) -mH_A{z). 

j&M " 


(32) 


We now have to examine Ylz'&E arbitrary partition A of S. To this end, we 

use ([71) and normalisation, followed by (I32p and a change of summation involving (|^ to calculate 


^ 7 ( 4 ( 2 :') = ^ A(z; y, x) [ 7 ( 4(2 + 4 - 4) “ 

^(iVr^ E ^B {RBi^))ix)4y)[HAi^ + ^x-Sy)-H^{z)] 


z'&E 


B& 0 ^ 2 (S) x,yeX 


E '’B 

i, _/Y' ^ ^ 

Be0^2{s) 

j£M 

E ''B 

E ® ^BU, - ^^)(^) 

B€0^2{S) 

j&M " 

E E 

E '’B ® Rj- 

isM j€0^2(Aj) 


B\.=J 

E E 

77 ® % - 7 ( 4 ) ( 2 ), 


R 


«b 


{N-m)\ - 

m H 


m 


A 


(z) 


j£M J7’€0^2(^j) 
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which agrees with (|3ip and proves the claim. 


□ 


We can now harvest some interesting consequences. First, Eq. m contains a meaningful 
expression for the derivative: 

Corollary 1. For A G ]P(5'), z & E, and the population process {Zt}t^o, we have 

^e^Ihaz,)\z» = z]\„„=Y. E ’>■(<■ n 

jeM j€0^2{Aj) 


The right-hand side has a plausible explanation. Namely, when block Aj splits into J ^ the other 
blocks in A retain their current type distribution (namely, (vr^^^ .z)). Independently of 

this, the parts of J pick their types from all individuals (with replacement), including those 
individuals that already carry other parts of AM\j, which is expressed by the tensor product 
with Rj{'k^,. z). 

Next, via ([26]) together with the fact that H_^{z) = [H_^{Zt) \ Zt = z], Eqns. ([271) and ([281) 

give rise to a system of differential equations for the expectations, namely: 

Corollary 2. For A G P(5) and the population process {Zt}t^Q, one has 

^E^[FF4(Z0]= j; 0^BE^[iFB(Zi)]. □ 

BeP(S) 

This will be the basis for our concrete calculations in the next section. 


7 Applications and examples 

Let us now apply our results to the cases of n = 2 and n = 3 sites. Expectations will always be 
with respect to ip, so we will abbreviate by E throughout. We will assume that Zq = z, i.e., 
that the initial population is deterministic. 


Two sites. Eor U = S = {1,2}, with the abbreviation r := the ODE system of 

Corollary [2] reads 

where we have dropped the upper index, which is always U. It follows that 

- %i},{2}})(^t)] = - E[(%i, 2}} -%i}.{2}})(^t)]- (34) 


Since L 


{{ 1 , 2 }} 


N-l 

N 


(FFj{i 2 }} — -^{{i},{ 2 }}), it follows that the expected two-point LDE also 
decays at ra te 2/N + r(N — 1)/N. In the c ase of two alleles per site, an equivalent formula has 
appeared in Bobrowski and Kimmell (120031 . Ex. 1). The correspo nding result in the diffusion 
limit goes back to lOhta and Kimural ( 1960l l. see also IPurrettl ( 20081 . Chap. 8.2). As noted there, 
it may seem surprising that the correlations also decay via resampling (even if r = 0); but recall 
that our Moran model with recombination is an absorbing Markov chain where a single type 
goes to fixation in the long run, that is, Zt will ultimately end up in a point measure. 
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The expected type distribution is now easily obtained from (1331) and (1341) via 

N-1 /■*. 




N 


E[(^^{{1,2}} - H{{l},{2}})iZr)]dT 


r(iV-l) 


N r{N -l) + 2 


1 — exp ( — 


r{N-l) + 2 


N 


t))E[(iJ|{i 2}} - i^{{i}_{2}})(^o)], 


where we have used that E 2 }}(-^o)] = ZqIN. The asymptotic behaviour is given by 

2 


lim E 

t^OO 


N 


_^0 I r(iV-l) 

2 + r(iV-l) iV 2 + r(iV-l) 


(Zo). 


(35) 


Since Zt will ultimately absorb in a point measure, we also know that 


lim E 

t—^oo 


N 


P[Zt absorbs in x] 5x, 




and thus P[Zt absorbs in x] = limt_,.oo P\Zt/N]{x) for all x € X. We can therefore read off the 
fixation probabilities from (1351) . With probability 2 +r{N-i) relative intensity of resampling), 

the type that wins is drawn from the initial distribution. With probability 2 +r{N^i) relative 
intensity of recombination), it is drawn from the distribution that results when the leading and 
the trailing segments are sampled from the initial population without replacement. 


Thre© sites. Now, consider U = S = {1, 2, 3}, together with the abbreviations ri = ^ 111112 , 3 }} 
and r 2 = 2 }{ 3 }}- Let us order the partitions of P(C/) as follows: 


{{1,2,3}} {{1},{2,3}} {{1,2},{3}} {{1,3},{2}} {{!},{2},{3}}. 

The generator of the partitioning process then reads 


/ -^iri+r2) ^r2 0 0 \ 




2 N-1^ 

N iV ^2 

2 

N jv^ *^2 

^r2 

iV2 

N-1 ^ 

4v2-^2 

(iV-l)(JV-2) 

iv2 *^2 


0 = 


N iV ' 1 


N n2 ' 1 


(iV-l)(Ar 2) 

• (36) 


2 

N 




1 L_7^(»-1+^2) 




\ 

0 

2 

N 

2 

N 

2 

N 

-- / 

N ' 



Recall that, by Corollary[51 we have ^ P[H{Zt)] = 0'E[H{Zt)], where H{Zt) := (■^t))^GP(U')- 

We now transform this system into a system in terms of correlation functions. Therefore, let 
L{Zt) = (L^(Xt))_ 4 gp(f/p From (|23|) . we know that L{Zt) = TH{Zt), where the transformation 
matrix is given by 



/ 1 

-1 

-1 

-1 

2 \ 


iV-2 


-1 

N-2 

-1 

N-2 


( iV - l )( Af - 2 ) 

W 

1 

iV-2 

1 

N-2 


1 

N-2 

-1 

1 

N-2 

1 

N-2 

1 

N-2 

N-2 

-1 


V 1 

1 

1 


1 / 


\ (iV-l)(iV-2) 

iV-2 

N-2 

N-2 


Consequently, ■^'E[L{Zt)] =T0T ^E[L(Zt)], where 


TQT-^ = 


N ' 

_ 2 _ 

N ' 



0 


0 


-E-E^ro 

N iV ' 2 


0 

L^(r-l+r2) 

0 

2 

N ' 

iV-1, 

N 

^ jv2 ^ hi+^ 2 ) 

0 


0 

■^(ri+r2) 

E_±ro 

N iV ' 2 

2 

N 
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0 o\ 

0 0 

ri 0 0 

-w-^Ci+^2) 0 
w-wCi+’'2) oj 


(37) 































In contrast to ()36p . the matrix TQT~^ has a nice subtriangular structure, from which we can 
already read off that the expected three-point LDE E[L||^ 2 , 3 }}(-^i)] (cf- (EH)) decays exponentially 
according to 


dt 




6A^+(iV-l)(iV-2)(ri +r 2 ) 

iV2 


E [-^{{1,2,3}} • 


As in the case of two sites, the decay rate contains contributions from resampling as well as 
from recombination. To extract more information, we recast TQT~^ into the diagonal form 
VTQT~^V~^ = D, where the entries of the diagonal matrix D are those on the diagonal of 
TQT~^, i.e., its eigenvalues. Consequently, 'E[L{Zt)] = DV~^ 'E[L{Zt)]. 

With the help of the subtriangular structure of T0T~^, the matrix V~^ can be calculated 
explicitly. It is again subtriangular, but somewhat unwieldy. To streamline the results, we now 
turn to the diffusion limit, with generator 0" of Definition El Then T and T~^ converge to 
matrices T" and (T")“^, respectively, with elements = p,{B, A) (T")^g = Sb^a^ 

A^B € ¥{U) (the latter is due to inverison from below). This yields 


T"0"{T")-^ 


( —(6+ei+e2) 0 0 0 0\ 

2 -(2+52) 0 0 0 \ 

2 0 -(2+5i) 0 0 , 

2 0 0 -(2+51+52) 0 / 

0 2 2 2 0 ^ 


where Qi = lim^v-ioo * = 1,2. Note that the rescaling of time has already been absorbed in 
0". In place of V~^, we now get 


(W')-' = 


1 

_ 2 _ 

(2+e2){4+ei) 

2 


(2+ei){4+e2) 

1 

2(2+ei+e2) 

4(ei e2+(2+gi +e2)(6+gi +g2)) 


0 

1 

2+^2 

0 

0 


0 

0 

1 

2+^1 

0 


°\ 

0 


^+Ql+Q2 

\ 4(gi g2 + (2+gI+g2)(6+gi+g2)) 2 2 2 -i , 

\ (2+e2)(2+£>2)(2+gi+g2)(®+^l + ^2) 2+^2 ’^+Q\-\-Q2 ' 


which diagonalises T0"T This shows that, in contrast to \U\ = 2, the linear combinations 
of E[L_ 4 (Zt)]’s that decay exponentially have coefficients depending on the recombination rates 
(with exception of E[L||i 2 , 3 }}(-^*)])- For example (4 +^> 1 ) E[L||i|| 2 , 3 }}(^*)] + 2E[L{|i 2 , 3 }}(^t)] 
is one such combination and decays at rate 2 -\- q^. Solution of the complete system is still 
possible due to the triangular structure; however, it is somewhat tedious since it involves the 
linear combination given in the last line of . Further progress may be possibl e if alternative 


scalin gs are employed, such as the loose linkage approach suggested recently by iJenkins et al 

(1201,^ 1. 


8 Conclusion 


Let us summarise our findings. We have described a marginal ancestral recombination process 
(ARP) and proved a duality result that relates the ARP with the Moran model forward in 
time, via so-called sampling functions. This was achieved by extending the recombinator 
formalism, which had previously proved useful in the context of deterministic recombination 
equations, to the stochastic setting. The ARP, together with th e duality result, reveals the 
genealogical structure hidden in the work of Bobrowski et al. ( 20101 ). who approached the matter 
by functional-analytic means and forward in time. It also leads to an explicit and closed system 
of ordinary differential equations for the expected sampling functions, from which the expected 
linkage disequilibria of all orders can be calculated. It is quite remarkable that such a closed ODE 
system exists: after all, the sampling functions are nonlinear, and the attempt to write down the 
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differential equation for the expect ation of a nonlinear quanti ty usually results in a hierarchy of 
equations that does not close; see iBaake and Hustedtl (j201l|) for more on the moment closure 
problem in the case of recombination. We would like to emphasise that the favourable structure 
is due to the marginalisation, which gives efficient access to correlation functions, but not to 
variances, for example. 

Unlike Bobrowski et al. ( 201Cll b we have not included mutation so far. However, since mutation 
acts independently of recombination, it should be straightforward to superimpose it on the 
population process as well as the partitioning process. It will be rewarding to study the interplay 
of mutation (which increases LDE) with recombination and resampling (which decrease LDE) 
within the framework established here. 
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